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Abstract 

We encode the binomials belonging to the toric ideal I a associated with an integral 
d x n matrix A using a short sum of rational functions as introduced by Barvinok 
Barvinok (1994); Barvinok and Woods (2003). Under the assumption that d,n are 
fixed, this representation allows us to compute the Graver basis and the reduced 
Grobner basis of the ideal I a, with respect to any term order, in time polynomial 
in the size of the input. We also derive a polynomial time algorithm for normal 
form computation which replaces in this new encoding the usual reductions typical 
of the division algorithm. We describe other applications, such as the computation 
of Hilbert series of normal semigroup rings, and we indicate further connections to 
integer programming and statistics. 

Key words: Grobner basis, toric ideals, Hilbert series, short rational function, 
Barvinok's algorithm, Ehrhart polynomial, lattice points, magic cubes and squares. 



1 Introduction 



In this note we present polynomial-time algorithms for computing with toric ideals 
and semigroup rings. For back ground on thes e algebraic objects a nd their interpla y 
with polyhedral geometry see ( Stanley , 1996f ) , ( Sturmfelsl . 1995 ) , ( Villarreal l200lh . 
Our results are a direct application of recent results by Barvinok and Woods (2003) 
on short encodings of rational generating functions (such as Hilbert series). 



Let A = ( ( 
hedron P 



'i j ■ 



be an integral d x n-matrix and 6 6 Z d such that the convex poly- 



{u e 



A ■ u = b and u > } is bounded. Barvinok (1994) 
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gave an algorithm for counting the lattice points in P in polynomial time when 
n — d is a constant. The input for Barvinok's algorithm is the binary encoding of 
the integers and fej, and the output is a formula for the multivariate generating 
function f(P) = J2a£Pnz n x<1 where x a is an abbreviation of This lone; 

polynomial with exponentially many monomials is encoded as a much shorter sum 
of rational functions of the form 

f(P) = Y±- r-, ; -. (1) 



Barvinok and Woods (2003) developed a set of powerful manipulation rules for us- 
ing these short rational functions in Boolean constructions on various sets of lattice 
points. In this note we apply their techniques to problems in combinatorial commu- 
tative algebra. Our first theorem concerns the computation of the toric ideal I a of 
the matrix A. This ideal is generated by all binomials x u — x v such that Au = Av. 
In what follows, we encode any set of binomials x u — x v in n variables as the formal 
sum of the corresponding monomials x u y v in 2n variables x±, . . . , x n , yi, . . . , y n . 

Theorem 1 Let A G Z dxn . Assuming that n and d are fixed, there is a polynomial 
time algorithm to compute a short rational function G which represents the reduced 
Grobner basis of the toric ideal I a with respect to any given term order -<. Given G 
and any monomial x a , the following tasks can be performed in polynomial time: 

(1) Decide whether x a is in normal form with respect to G. 

(2) Perform one step of the division algorithm modulo G. 

(3) Compute the normal form of x a modulo the Grobner basis G. 

Our research group at UC Davis has developed a computer program, called LattE, 
which efficiently counts the lattice points in any rational polytope by computing 
its Barvinok representation (1). The Grobner basis and normal form algorithms of 
Theorem 1 are currently being implemented in LattE. It is important to note that 
the Grobner basis G which will be output by LattE is a rational function. It 
is not the long list of binomials produced by all other computer algebra systems. 

Example 2 Fix n = 4, d = 2 and let m > 3 be an arbitrary integer. Consider 

m m — 1 1 

inputting the matrix A = and the lexicographic term order 

1 m — 1 m 
into LattE. The task is to compute the kernel I a of 

k[xi,x 2 , x 3 , X4] — > k[s, t] , X\ 1 — ^ s m , x 2 1— > s m ~ 1 t, x 3 1— > st m ~ l , x 4 1— » t m . 

The output produced by LattE would consist of the rational function 

2(l \m—l 1 \m— 1 



G(x,y) = x x x A y 2 yz + x 2 x™ y™ + 



xi x 3 y 2 " [ [x 1 y 2 ) - [x 3 y A ) 



xi y 2 - x 3 1/4 
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This rational function is a polynomial whose number of terms is m + 1 and hence 
grows exponentially in the size of the input. Yet, the running time for computing 
G(x, y) is bounded by a polynomial in log(m). It is an interesting exercise to perform 
the tasks (1), (2) and (3) in Theorem 1 for G(x,y) and the monomial 

The proof of Theorem 1 will be given in Section 2. Sp ecial attention will be paid to 
the Projection Theorem ( Barvinok and Woodsl 20031 Theorem 1.7) since projection 



of short rational functions is the most difficult step to implement. Its practical 
efficiency has yet to be investigated. Our proof of Theorem 1 does use the Projection 
Theorem, but our Proposition 8 in Section 2 shows that a non-reduced Grobner basis 
can be computed in polynomial time without using the Projection Theorem. 

In Section 3 we present wha t we call the homogeni zed Barvinok algorithm. This algo- 
rithm was first outlined in ( De Loera et al. . 2003) and it was recently implemented 



in LattE. Like the original version in ( BarvinokT 1994 ) . it runs in polynomial time 



when the dimension is fixed. But it performs much better in practice (1) when com- 
puting the Ehrhart series of polytopes with few facets but many vertices; (2) when 
computing the Hilbert series of normal semigroup rings. We show its effectiveness 
by solving the classical counting problems for 5 x 5 magic squares (all row, column 
and diagonal sums are equal) and 3x3x3x3 magic hypercubes (Theorem 12). 

A normal semigroup S is the intersection of the lattice Z n with a rational convex 
polyhedral cone in K n . The Hilbert series of S is the rational generating function 
yj agS x a . Barvinok and Woods (2003) showed that this Hilbert series can be com- 
puted as a short rational generating function. We show that this computation can be 
done without the Projection Theorem when the semigroup is known to be normal. 

Theorem 3 Under the hypothesis that the ambient dimension n is fixed, 

1 ) the Ehrhart series of a rational convex polytope given by linear inequalities can be 
computed in polynomial time. The Projection Theorem is not used in the algorithm. 

2) The same applies to computing the Hilbert series of a normal semigroup S. 

In the final section of the paper we sketch applications of the above algebraic theory 
to Integer Programming and Statistics. These results will be explored in detailed 
elsewhere. 



2 Computing Toric Ideals 



We assum e that the reader is familia r with toric ideals and Grobner bases as pre- 
sented in (|Oox et all [l992t ISturmfeleJ . Il995l) . Barvinok and Woods (2003) showed: 



Lemma 4 (Theorem 3.6 in (|Barvinok and Woodsl 20031 )^ LetSi,Si be finite 
subsets of 7U 1 , for n fixed. Let f{S\,x) and f(S 2 ,x) be their generating functions, 
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given as short rational functions with at most k binomials in each denominator. 
Then there exist a polynomial time algorithm, which, given f(Si,x), computes 



f(Sif]S 2 ,x) = 7j • — r 

f£ ' (1-X v *)...(l-X v i 



with s < 2k, where the 7$ are rational numbers, and Ui,Vij nonzero integers. 

We will use this Intersection Lemma to extract special monomials present in the 
expansion of a generating function. The essential step in the int ersection algorithm 
is the use of the Hadamard product ( Barvinok and Woodsl 200.4 Definition 3.2) and 



a special monomial substitution. The Hadamard product is a bilinear operation on 
rational functions (we denote it by *). The computation is carried out for pairs of 
summands as in (1). Note that the Hadamard product mi * m 2 of two monomials 
m 1; m 2 is zero unless m\ = m 2 . We present an example of computing intersections. 

Example 5 Let Si = { x G R : % — 2 < x < % } D Z for % — 1, 2. We rewrite their ra- 
tiona l generating functions as in the proof of Theorem 3.6 in ( Barvinok and Woodsl 



200i): f(S u z) = + (Yzf^y = ^=^y + {izjpiy = 9u + 9u, and f(S 2 ,z 

1 z 2 —z^ 1 z 2 

(l=i) + (l- z -l) = (l- 2 -i) + = + 922- 

We need to compute four Hadamard products between rational functions whose de- 
no minators are products of bin omials and denominators are monomials. Lemma 3.4 
in 



Barvinok and Woodsl (l2003f ) says that, for our example, these Hadamard prod- 



ucts are essentially the same as computing the functions (1) of the auxiliary poly- 
hedron {(ei,e 2 )|pi + ai€\ = p 2 + 0,2^2, Q > 0} where pi,p 2 are the exponent of 
numerators of g'^s involved and a\, a 2 are the exponents of the binomial denomina- 
tors. For example, the Hadamard product gn * g 22 corresponds to the polyhedron 
{(ei, e 2 )|e 2 = 4 + ei, 6j > 0}. The contribution of this half line is — hz^=t\ ■ We find 

f(S 1 ,z)*f(S 2 ,z) = g n * g 2 i + g 12 * g 22 + g 12 * g 2 \ + gu * g 22 

z~ 2 z z~ x z~ 2 



z- 1 ) (l-z- 1 ) (l-z- 1 ) (l-z- 1 ) 

= i + z = ffans^z). 



z — z 1 



z- 1 



Another key subroutine introduced by Barvinok and Woods is the following Projec- 
tion Theorem. In both Lemmas 4 and 6, the dimension n is assumed to be fixed. 



Lemma 6 (Theorem 1.7 in ( Barvinok and Woodsl 20031 )^ Assume the dimen 



sion n is a fixed constant. Consider a rational polytope PcK" and a linear map 
T : Z n — > Z fc . There is a polynomial time algorithm which computes a short repre- 
sentation of the generating function f(T(P D 7* n ),x). 

We repre sent a term order -< on mon omials in x\, . . . , x n by an integral n x n-matrix 
W as in ( Mora and Robbianol fl998h . Two monomials satisfy x a -< x 13 if and only 
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if Wet is lexicographically smaller than Wf3. In other words, if w±, . . . , w n denote 
the rows of W, there is some j G {1, ... ,71} such that WiCt = Wi/3 for z < j, and 
WjQ; < Wjf3. For example, W = I n describes the lexicographic term ordering. In 
what follows, we will denote by -<w the term ordering defined by W. 

Lemma 7 Let S C Z" be finite. Suppose the polynomial f(S,x) = J2p£S xf3 
represented as a short rational function and let -<w be a term order. We can extract 
the (unique) leading monomial of f(S,x) with respect to -<w, i n polynomial time. 

Proof: The term order -<w is represented by an integer matrix W . For each of the 
rows Wj of W we perform a monomial substitution Xj := xjt™ 31 . Such a monomial 



substitution can be computed in polynomial time by (jBarvinok and Woodsl 20031 
Theorem 2.6). The effect is that the polynomial f(S, x) gets replaced by a polynomial 
in the t and the x's. After each substitution we determine the degree in t. This is 
done as follows: We want to do calculations in univariate polynomials since this is 
faster so we consider the polynomial g(t) = f(S, l,t), where all variables except t 
are set to the constant one. Clearly the degree of g(t) in t is the same as the degree 
of f(S,x',t). We create the interval polynomial i\p, q ]{t) = J2l= p t l which obviously 
has a short rational function representation. Compute the Hadamard product of 
and ifp q ] with g(t). This yields those monomials whose degree in the variable t lies 
between p and q. We will keep shrinking the interval [p, q] until we find the degree. 
We need a bound for the degree in t of g(t) to start a binary searc h. A polynomia l 



upper bound U can be found via the estimate in Theorem 3.1 of ( LasserreL l2003f ) 



by easy manipulation of the numerator and denominator of the fractions in g(t). In 
no more than log{U) steps one can determine the degree in t of f(S,x,t) by using 
a standard binary search algorithm. 

Once the degree r in t is known, we compute the Hadamard product of f(S,x,t) 
and i[ r>r ], and then compute the limit as t approaches 1. This can be done in poly- 
nomial time using residue techniques. The limit represents the subseries H(S, x) = 
Y,0- W] = r x ^ ■ Repeat the monomial and highest degree search for the row Wj + i,Wj + 2, 
etc. Since -<w is a term order, after doing this n times we will have only one single 
monomial left, the desired leading monomial. □ 

Proposition 8 Let A E Z dxn , W G Z nxn specifying a term order -< w , and assume 
that d and n are fixed. 

1 ) There is a polynomial time algorithm to compute a short rational function G 
which represents a universal Grobner basis of 1a- 

2) Given the term order -<\y and a short rational function encoding a (possibly infi- 
nite) set of binomials ^x u y v , one can compute in polynomial time a short rational 
function encoding only those binomials such that x v -<w x u . 

3) Suppose we are given a sum of short rational functions fix) which is identical, 
in a monomial expansion, to a single monomial x a . Then in polynomial time we can 
recover the (unique) exponent vector a. 
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Proof: 1) Denote by Wi the i-th row of the matrix W which specifies the term order. 
Set M — (d+l)(n — d)D(A) where D(A) is the large st absolut e value of any dx d- 



sub determinant of A. Using Barvinok's algorithm in (jBarvinokl . 119941 ) . we compute 
the following generating function in In variables: 

G(x,y) = x u y v : Au = Av and < Ui, ih < M}. 

This is the sum over all lattice points in a rational polytope. Lemma 4.1 and Theorem 
4.7 in Chapter 4 of (ISturmfelsl . ll995h imply that the toric ideal I a is generated by the 



finite set of binomials x u — x v corresponding to the terms x u y v in G(x, y). Moreover, 
these binomials are a universal Grobner basis of Ia- 

2) Suppose we are given a short rational generating function G (x,y) = J2 xU y v 
representing a set of binomials x u — x v in I a, for instance Go = G in part (1). In the 
following steps, we will alter the series so that a term x u y v gets removed whenever 
u is not bigger than v in the term order -<w Starting with H = Gq, we perform 
Hadamard products with short rational functions f(S; x, y) for S C Z 2n . 

Set Hi = f({(u,v) : w(a = Wiv}), and G { = H i ^ 1 *f({(u,v) : w^ > WiV + 1}). 

All monomials x u y v G Gj have the property that w{ii = WiV for i < j, wju > WjV, 
and thus v -<w u - On the other hand, if v -<w u then there is some j such that 
WiU = WiV for i < j, WjU > WjV, and we can conclude that x u y v G Gj. This 
proves that G = G\ U G<i U . . . U G n encodes exactly those binomials in Gq that are 
correctly ordered with respect to -<w We have proved our claim since all of the 
above constructions can be done in polynomial time. 

3) Given f(x) we can compute in polynomial time the partial derivative df(x)/dxi. 
This puts the exponent of Xi as a coefficient of the unique monomial. To compute the 
derivative can be done in polynomial time by the quotient and product derivative 
rules. Each time we differentiate a short rational function of the form 

x bi 



(1 - £ c m)(1 - x C2 > ! ) ... (1 - x Cd >>) 

we add polynomially many (binomial) factors to the numerator. The factors in the 
numerators should be expanded into monomials to have again summands in short 
rational canonical form 7- — jt- w- — T . — ettt- Note that at most 2 d monomials 

(1— x L ' l )(l — x ..(1— x <M) 

appear (a constant) each time. Finally, if we take the limit when all variables Xi go 
to one we will get the desired exponent. □ 

Example 9 Using LattE we compute the set of all binomials of degree less than or 

1111 

equal 10000 in the toric ideal I a of the matrix A = . This matrix repre- 

12 3 

sents the Twisted Cubic Curve in algebraic geometry. We find that there are exactly 
195281738790588958143425 such binomials. Each binomial is encoded as a monomial 
x^x^x^x^y^y^yTvT ■ The computation takes about 40 seconds. The output is a 
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sum of 538 simple rational functions of the form a monomial divided by a product 
such as (l - sg) (l - (1 - xiyi) (1 - x 1 x 3 y 2 2 ) (1 - x 3 y 3 ) (1 - x 2 y 2 ). □ 

Proof of Theorem 1: Proposition 8 gives a Grobner basis for the toric ideal I a in 
polynomial time. We now show how to get the reduced Grobner basis. 

Step 1. Compute the generating function which encodes all binomials in I a- 

f(x,y) = ^2{x u y v : Au = Av and u, v > }, 

This computation is similar to part 1 of Proposition 8 except that there is no upper 
bound M. Next we wish to remove from f(x,y) all incorrectly ordered binomials 
(i.e. those monomials x u y v with u -<w v instead of the other way around). We do 
this following part 2 of Proposition 8. Abusing notation let us still call f(x,y) the 
resulting sum of rational functions. Let now g(x) be the projection of f(x,y) onto 
the first group of variables. Thus g(x) is the sum over all non-standard monomials, 
and it can be computed in polynomial time by Lemma 6. 

Step 2. Write = n f° r the generating function of all x- monomials. We 

i=i % 

compute the following Hadamard product of n rational functions in x: 

(r^r ~ Xl ' * (r^r ~ X2 ' 9<yX ^) * " * * (t^x ~ Xn ' 

This is the generating function over those monomials all of whose proper factors are 
standard monomials modulo the toric ideal I a- 

Step 3. Let h(x, y) denote the ordinary product of the result of Step 2 with 

g(y) = ~^2{y v ■ v standard monomial modulo I a }• 

Thus h(x, y) is the sum of all monomials x u y v such that x v is standard and x u 
is minimally non-standard. Compute the Hadamard product G(x,y) := f(x,y) * 
h(x, y). This is a short rational representation of a polynomial, namely, it is the sum 
over all monomials x u y v such that the binomial x u — x v is in the reduced Grobner 
basis of I a with respect to W. 

We next prove claims 1 and 2. Let G(x, y) be the reduced Grobner basis of I a encoded 
by the rational function above, and let M be the degree bound of Proposition 8. Let 
b(x,y) be the rational function representing {(u, v) : < u < a, < v < M}. The 
Hadamard product G(x, y) = G(x, y)*b(x, y) is computable in polynomial time and 
encodes exactly those binomials in G that can reduce x a . If G is empty then x a is in 
normal form already, otherwise we use Lemma 7 to find an element (u, v) £ G and 
reduce x a to x a ~ u+v . 

It is worth noting that analytic calculations may now be used as part of algebraic 
algorithms: Suppose again we wish to decide whether x a is in reduced normal form 
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or not. Take G(x,y) as before and compute F(x) = G(x, 1). This can be done using 
monomial substitution ( Barvinok and Woods! . 2003f ). Next compute the integral 



)»/"'/ 



x^---x- a -)F(x) ^ 



(2ttz)™ / / (1- Xl )...(l- Xn 

Here < e±, . . . , < 1 are different numbers such that we can expand all the jz^ 
into the power series about 0. It is possible to do a partial fraction decomposition of 
the integrand into a sum of simple fractions. The integral is a non-negative integer: 
it is the number of ways that the monomial x a can be written in terms of the leading 
monomials of the Grobner bases G. 

We now present the algorithm for claim 3 in Theorem 1. A curious byproduct of 
representing Grobner bases with short rational functions is that the reduction to 
normal form need not be done by dividing several times anymore: 

Step 4. Let f(x,y) and g(x) as above and compute the Hadamard product 

H(x,y) := f( x ,y)J(J—).(J— g( y ))). 

\1 — x 1 — y J 

This is the sum over all monomials x u y v where x v is the normal form of x u . 

Step 5. We use H(x, y) as one would use a traditional Grobner basis of the ideal I a- 
The normal form of a monomial x a is computed by forming the Hadamard product 

H(x,y) * 



Since this is strictly speaking a sum of rational functions equal to a single monomial, 
applying Part 3 of Proposition 8 concludes the proof of Theorem 1. □ 



3 Computing Normal Semigroup Rings 



We observed in ( De Loera et al 



th at a major practical bottleneck of the 



original Barvinok algorithm in ([Barvinokl . [1994J) is the fact that a polytope may 
have too many vertices. Since originally one visits each vertex to compute a rational 
function at each tangent cone, the result can be costly. For example, the well-known 
polytope of semi-magic cubes in the 4x4x4 case has over two million vertices, 
but only 64 linear inequalities describe the polytope. In such cases we propose a 
homogenization of Barvinok's algorithm working with a single cone. 

There is a second motivation for looking at the homogenization. Barvinok and Woods 
( Barvinok and Woodsl l2003h proved that the Hilbert series of semigroup rings can 
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be computed in polynomial time. We show that for normal semigroup rings this can 
be done simpler and more directly, without using the Projection Theorem. 



Given a rational polytope P in M. n 1 , we set i(P,m) = #{z £ Z n 1 : z £ mP}. The 
Ehrhart series of P is the generating function Y^m=o rn)t m . 

Input: A full-dimensional, rational convex polytope P in M n_1 specified by linear 
inequalities and linear equations. 
Output: The Ehrhart series of P. 

(1) Place the polytope P into the hyperplane defined by x n = 1 in IR n . Let K be 
the n-dimensional cone over P, that is, K = cone({(p, 1) : p £ P}). 

(2) Compute the polar cone if*. The normal vectors of the facets of K are exactly 
the extreme rays of K*. If the polytope P has far fewer facets then vertices, 
then the number of rays of the cone K* is small. 

(3) Apply Barvinok's decomposition of K* into unimodular cones. Polar- 
i ze back each of these cones. It is known, e.g. Corollary 2.8 in 
(jfjarvinok and Pommersheiml . 19991 ). that by dualizing back we get a unimod- 



ular cone decomposition of K. All these cones have the same dimension as K. 
Retrieve a signed sum of multivariate rational functions which represents the 
series Daeifnz™ x °" ■ 

(4) Replace the variables xt by 1 for % < n — 1 and output the resulting series in 
t — x n . This can be done using the methods in ( De Loera et al. . 2003j ). 



We recall that one of the key steps in Barvinok's algorithm is that any cone can be 
decomposed as the signed sum of (indicator functions of) unimodular cones. 



Theorem 10 (see (jBarvinokL Il994h ^ Fix the dimension n. Then there exists a 



polynomial time algorithm which decomposes a rational polyhedral cone K C W 1 into 
unimodular cones Ki with numbers £ { — 1, 1} such that 

f(K) = J2af( K d, \I\<oo. 



Originally, Barvinok had pasted together such formulas, one for each vertex of a 
polytope, using a result of Brion. The point is that this can be avoided: 

Proof of Theorem 3: We first prove part (1). The algorithm solving the problems 
is Algorithm 3. Steps 1 and 2 are polynomial when the dimension is fixed. Step 3 
follows from Theorem 10. We require a special monom ial substitution, possibly wit h 
some poles. This can be done in polynomial time by ( Barvinok and Woodsl 2003| ). 



Part (2): From the characterization of the integral closure of the semigroup S as the 
intersection of a pointed polyhedral cone with the lattice Z™ is clear that Algorithm 
1, with the modification that the cone K in question is given by the rays of the cone 
(the generators of the monomial algebra). In fixed dimension one can transfer from 
the the extreme rays representation of the cone or to the facet representation of the 
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cone in polynomial time. □ 



Corollary 11 Given a normal semigroup ring R of fixed Krull dimension, there is 
a polynomial time algorithm which decides whether R is Gorenstein. 

Proof: Let R be a normal semigroup ring for a semigroup of Z n . Hoc hster's theorem 
says that the normal semigroup rings are Cohen-Macaulay domains ( Stanley . 1996| ). 



Denote by F(R, z) the generating function of the monomials of normal semigroup 
ring R (computable in polynomial time by previous theorem). Then by Theorem 
12.7 in rtStanlevLllQQfih . it is enough to check that F(R,z) = (-l) n z a F(R, l/z) for 
some a G Z n efficiently The change of variables can be done in polynomial time and 
thus get F(R, l/z). To check whether F(R,z)/F(R, l/z) is a single polynomial we 
can compute the monomial evaluation z% = 1 for % — 1 . . . d. □ 

Each pointed afhne semigroup S C Z™ can be graded. This means that there is a 
linear map deg : S — > N with deg(x) = if and only if x = 0. Given a pointed graded 
afhne semigroup define S r to be the set of all degree r elements, i.e. S r = {x G S : 
deg(x) = r}. The Hilbert series of S is the formal power series Hs(t) = J2kLo #(Sr)t r - 
Algebraically, this is just the Hilbert series of the semigroup ring C[S]. It is a well- 
known property that H$ is represented by a rational function of the form 

Q(t) 



(i-t*)(i-f da ).-.(i-t dn 



wher e Q(t) is a polynomial of degree less than d\ + . . . + d n (see Chapter 4 (jStanlev 



Il997t)). Several other met hods had been tried to compute the Hilbert series explicitly 



see (jAhmed et al. . 2003| ) for references). One of the most well-known challenges was 



that of counting the 5x5 magic squares of magic sum n. Similarly several authors 
had tried before to compute the Hilbert series of the 3x3x3x3 semi-magic 
cubes. It is not difficult to see this is equivalent to determining an Ehrhart series. 
Using Algorithm 1 we finally present the solution, which had been inaccessible using 
Grobner bases methods. For comparison, the reader familiar with Grobner bases 
computations should be aware that the 5x5 magic squares problem required a 
computation of a Grobner bases of a toric ideal of a matrix A with 25 rows and over 
4828 columns. Our attempts to handle this problem with CoCoA and Macaulay2 were 
unsuccessful. We now give the numerator and then the denominator of the rational 
functions computed with the software LattE: 

Theorem 12 

The generating function J2n>o f( n )t n f or the number fin) of 5 x 5 magic squares 
of magic sum n is given by the rational function p(t) / q(t) with denominator 

p(t) = t 76 +28t 75 +639 t 74 +11050t 73 + 136266t 72 + 125 5 8 33 t 71 +9120009t 70 +54389347t 69 + 
274778754 1 68 + 1204206107 1 67 + 4663304831 1 66 + 16193751710 1 65 + 51030919095 1 64 + 
147368813970 t 63 + 393197605792 t 62 + 975980866856 t 61 + 2266977091533 t 60 + 
4952467350549 1 59 + 10220353765317 1 58 + 20000425620982 1 57 + 37238997469701 1 56 + 
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66164771134709 1 55 + 112476891429452 1 54 + 183365550921732 t 53 + 287269293973236 t b2 + 
433289919534912 1 51 +630230390692834 1 50 +885291593024017 1 49 +1202550133880678 t 48 + 
1581424159799051 t 47 + 2015395674628040 t 46 + 2491275358809867 t 45 + 
2989255690350053 i 44 + 3483898479782320 t 43 + 3946056312532923 t 42 + 
4345559454316341 t 41 + 4654344257066635 t 40 + 4849590327731195 t 39 + 
4916398325176454 t 38 + 4849590327731195 t 37 + 4654344257066635 t m + 
4345559454316341 t 35 + 3946056312532923 t 34 + 3483898479782320 t 33 + 
2989255690350053 t 32 + 2491275358809867 t 31 + 2015395674628040 t 30 + 
1581424159799051 t 29 + 1202550133880678 t 28 + 885291593024017 t 27 + 
630230390692834 1 26 +433289919534912 1 25 +287269293973236 1 24 + 183365550921732 1 23 + 
112476891429452 t 22 + 66164771134709 t 21 + 37238997469701 t 20 + 20000425620982 t 19 + 
10220353765317 t 18 + 4952467350549 t 17 + 2266977091533 t w + 975980866856 t 15 + 
393197605792 1 14 + 147368813970 1 13 +51030919095 1 12 +16193751710 1 11 +4663304831 1 10 + 
1204206107 1 9 + 274778754 1 8 + 54389347 1 7 + 9120009 1 6 + 1255833 1 5 + 136266 1 4 + 11050 1 3 + 
639 1 2 + 28 1 + 1 and numerator 

q (t) = {t 2 - l) 10 {t 2 + t + l) 7 (t 7 - l) 2 (t 6 + t 3 + 1) (t 5 + t 3 + t 2 + t + l) 4 (1 - t) 3 {t 2 + l) 4 . 

T7ie generating function J2 n >o fi n )t n f or ^ e number f(n) of 3x3x3x3 magic 
cubes with magic sum n is given the rational function r(t)/s(t) where 

i 54 + 150 1 51 + 5837 1 48 + 63127 1 45 + 331124 1 42 + 1056374 1 39 + 2326380 1 36 + 3842273 1 33 + 
5055138 i 30 +5512456 1 27 +5055138 1 24 +3842273 1 21 +2326380 i 18 +1056374 1 15 +331124 1 12 + 
63127 1 9 + 5837 1 6 + 150 1 3 + 1 and 

q (t) = (t 3 + l) 4 (t 12 + t 9 + t 6 + t 3 + 1) (1 - t 3 ) 9 (t 6 + t 3 + 1) . 



4 Applications 



As explained in Chapter 5 of the book Sturmfelsl ( 1995| ). Grobner bases can be 



useful in the context of integer programming, serving as universal test sets of fam- 
ilies of integer programs, and in statistics, where they can be used as the Markov 
basis moves used to generated elements uniformly at random (e.g contingency ta- 
bles counting). Therefore the fact that we can compute Grobner bases and normal 
forms in polynomial time (under certain hypothesis) can then be used to prove the 
following results: 

Corollary 13 Let A e Z dxn , b G Z d , W e Z n . Assume that d and n are fixed. 
There is a polynomial time algorithm to solve the integer programming problem 
min X (zp n znWx where P(b) = {x\Ax = b, x > 0}. 

sketch of proof: Make the cost vector W into a term order by breaking ties of 
the order mi > mp. if Wmi > W m^,. You can do this via lexicographic ordering. 
From Chapter 5 of Sturmfelsl ( 19951 ) the integral optimum of P can be obtained from 



the Grobner basis obtained in Theorem 1 and then computing the normal form of 
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Table 1 

The conditions for retinoblastoma RB1-VNTR genotype data from the Ceph database. 

the monomial x b with respect to the Grobner basis. Since both steps can be done 
efficiently the corollary follows. 

Another application is to the uniform sampling of lattice points inside polyhedra 
of the form P(b) = {x G M^Ae = b, x > 0}. Given M be a finite set such that 
M C {x G 7i d \Ax = 0}. We define the graph Gb such that its nodes are all the lattice 
points inside of P and there is an undirected edge between a node u and a node v 
iff it — v G M. In general this graph may not be connected for arbitrary choices of 
M. We say M is a Markov basis if Gb is a connected graph for all b. 

Corollary 14 Given A G Z dxn , where d and n are fixed, there is a polynomial time 
algorithm to compute a multivariate rational generating function for a Markov basis 
M associated to A. This is presented as a short sum of rational functions. 

We conclude with another with numeric question. Ian Dinwoodie communicated to 
us the problem of counting all 7 x 7 contingency tables whose entries are nonnegative 
integers Xi, with diagonal entries multiplied by a constant as presented in Table 1. 
The row sums and column sums of the entries are given there too. Using LattE we 
obtained the exact answer 8813835312287964978891 
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